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Abstract 

In this paper wc address the problem of uncertainty management 
for robust design, and verification of large dynamic networks whose 
performance is affected by an equally large number of uncertain pa- 
rameters. Many such networks (e.g. power, thermal and communica- 
tion networks) are often composed of weakly interacting subnetworks. 
We propose intrusive and non-intrusive iterative schemes that exploit 
such weak interconnections to overcome dimensionality ciirsc associ- 
ated with traditional uncertainty quantification methods (e.g. gen- 
eralized Polynomial Chaos, Probabilistic Collocation) and accelerate 
uncertainty propagation in systems with large number of uncertain 
parameters. This approach relies on integrating graph theoretic meth- 
ods and waveform relaxation with generalized Polynomial Chaos, and 
Probabilistic Collocation, rendering these techniques scalable. Wc an- 
alyze convergence properties of this scheme and illustrate it on several 
examples. 

1 Introduction 

The issue of management of uncertainty for robust system operation is of 
interest in a large family of complex networked systems. Examples include 
power, thermal and communication networks which arise in several instances 
such as more electric aircrafts, integrated building systems and sensor net- 
works. Such systems typically involve a large number of heterogeneous, 
connected components, whose dynamics is affected by possibly an equally 
large number of uncertain parameters and disturbances. 

Uncertainty Quantification (UQ) methods provide means of calculating 
probability distribution of system outputs, given probability distribution of 



input parameters. Outputs of interest could include for example, latency 
in communication network, power quality and stability of power networks, 
and energy usage in thermal networks. The standard UQ methods such as 
Monte Carlo(MC) |lj either exhibit poor convergence rates or others such 
as Quasi Monte Carlo (QMC) [2] [3], generalized Polynomial Chaos(PC) ^ 
and the associated Probabilistic Collocation method (PCM) [5], suffer from 
the curse of dimensionality (in parameter space) , and become practically in- 
feasible when applied to network as a whole. Improving these techniques to 
alleviate the curse of dimensionality is an active area of current research (see 
[6] and references therein): notable methods include sparse-grid collocation 
method [Z],[S! and ANOVA decomposition \W\ for sensitivity analysis and 
dimensional reduction of the uncertain parametric space. However, none of 
such extension exploits the underlying structure and dynamics of the net- 
worked systems. In fact, many networks of interest (e.g. power, thermal and 
communication networks) , are often composed of weakly interacting subsys- 
tems. As a result, it is plausible to simplify and accelerate the simulation, 
analysis and uncertainty propagation in such systems by suitably decom- 
posing them. For example, authors in [lOl [H] used graph decomposition 
to facilitate stability and robustness analysis of large-scale interconnected 
dynamical systems. Mezic et al. [Ll2j used graph decomposition in con- 
junction with Perron Frobenius operator theory to simplify the invariant 
measure computation and uncertainty quantification, for a particular class 
of networks. While these approaches exploit the underlying structure of the 
system, they do not take advantage of the weakly coupled dynamics of the 
subsystems. 

In this paper, we propose an iterative UQ approach that exploits the 
weak interactions among subsystems in a networked system to overcome 
the dimensionality curse associated with traditional UQ methods. We refer 
to this approach as Probabilistic Waveform Relaxation (PWR), and propose 
both intrusive and non-intrusive forms of PWR. PWR relies on integrating 
graph decomposition techniques and waveform relaxation scheme, with gPC 
and PCM. Graph decomposition to identify weakly interacting subsystems, 
can be realized by spectral graph theoretic techniques jl3j.|14j. Waveform 
relaxation [15] (WR), a parallelizable iterative method, on the other hand, 
exploits this decomposition and evolves each subsystem forward in time in- 
dependently but coupled with the other subsystems through their solutions 
from the previous iteration. In the intrusive PWR, the subsystems obtained 
from decomposing the original system are used to impose a decomposition 
on system obtained by Galerkin projection based on the gPC expansion. 
Further the weak interactions are used to discard terms which are expected 
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to be insignificant in the gPC expansion, leading to what we call an Ap- 
proximate Galerkin Projected (AGP) system. We then propose to apply 
WR relaxation on the decomposed AGP system to accelerate the UQ com- 
putation. In the non-intrusive form of PWR, rather than deriving the AGP 
system, one works directly with subsystems obtained from decomposing the 
original system. At each waveform relaxation iteration we propose to apply 
PCM at subsystem level, and use gPC to propagate the uncertainty among 
the subsystems. Since UQ methods are applied to relatively simpler sub- 
systems which typically involve a few parameters, this renders a scalable 
non- intrusive iterative approach to UQ. We prove convergence of the PWR 
approach under very general conditions. Note that spectral graph decom- 
position can be done completely in a distributed fashion using a recently 
developed wave equation based clustering method Moreover, one can 
further exploit timescale separation in the system to accelerate WR using 
an adaptive form of WR [T7]. PWR when combined with wave equation 
based distributed clustering and adaptive WR can lead to highly scalable 
and computationally efficient approach to UQ in complex networks. 

This paper is organized in six sections. In section [2] we give the math- 
ematical preliminaries for setting up the UQ problem for networked dy- 
namical systems, and present an overview of gPC and PCM techniques. In 
section |3] we discuss graph decomposition and waveform relaxation methods, 
which form basic ingredients of PWR. Here we also describe adaptive WR 
and wave equation based distributed graph decomposition techniques. We 
introduce the intrusive and non-intrusive PWR in section [4] through a simple 
example, and then describe these methods in a more general setting. We 
also prove convergence of PWR, and analyze the scalability of the method. 
In section [5] we illustrate the intrusive and non-intrusive PWR on several 
examples. Finally, in section [6] we summarize the main results of this paper, 
and present some future research directions. 

2 Uncertainty Quantification in Networked Sys- 
tems 

Consider a nonlinear system described by a system of random differential 
equation 

xi = /i(x,^i,t), 

Xn = /n(x,^ri)i)) (1) 
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where, f = (/i, /2, • • • , fn) G M" is a smooth vector field, x = {xi,X2, • • • , Xn) G 
M" are state variables, G i?^' is vector of random variables affecting the 
i—th. system. Let ^ = {^f , ■ ■ ■ ,Cn)"^ ^ be the p = Yll=iPi dimensional 
random vector of uncertain parameters affecting the complete system. The 
solution to initial value problem x(to) = xq will be denoted by x(t; ^), where 
for brevity we have suppressed the dependence of solution on initial time 
and initial condition xq. We shall assume that the system ([T]), is Lipschitz 

||f(xi,^,t)-f(x2,e,t)|| <L(e)||xi-X2||, (2) 

where, the Lipschitz constant -L(^) depends on the random parameter vector 
and II • II is a Euclidean norm. We will assume that supgg]gp L{C) = L < oo. 
Let us also define a set of quantities 

z = (.zi,Z2,-- ■ ,Zd) = G(x) = (51 (x), • • • ,5d(x)), (3) 

as observables or quantities of interests. The goal is to numerically establish 
the effect of input uncertainty of ^ on output observables z. Naturally, the 
solution for system ([T]) and the observables ([3| are functions of same set of 
random variables ^, i.e 

x = x(t;e), z = z(t,e) = G(x). (4) 

In what follows we will adopt a probabilistic framework and model 
^ = (Ci)C2> • • • j'^p) as a p— variate random vector with independent com- 
ponents in the probability space {Q,A,V), whose event space is Q and is 
equipped with a— algebra A and probability measure V. Throughout this 
paper, we will assume that the parameters S = {Ci,^2r"" jCp} are mu- 
tually independent of each other. Let pj : Fj — ?■ M"*" be the probability 
density of the random variable ^j(a;), with Fj = C M being its im- 

age. Then, the joint probability density of any random parameter subset 
A = {^h,^i2,- ■ ■ , Cim} C S is given by 

|A| 

PAi(.ii,--- = YlpijiCi,), V(^ii,--- ,^i,„) e La, (5) 
with a support 

|A| 

rA = nri. cmI^i, (6) 

where, | • | denotes the cardinality of the set. Without loss of generality we 
will assume that Fj = [—1 1], i = 1, • • • ,p. 
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Remark 2.1. While throughout the paper we will work with ODEs ^ with 
parametric uncertainty, the PWR framework developed in this paper can he 
naturally extended to deal with 1 ) System of differential algebraic equations 
(DAEs), and 2) Time varying uncertainty. Both these extensions would he 
illustrated through examples in section\^ 

2.1 Uncertainty Quantification Methods 

In this section, we describe two interrelated UQ approaches: generahzed 
polynomial chaos (gPC) and probabilistic collocation method (PCM). The 
gPC is an intrusive approach which requires explicit access to system equa- 
tions ([l]), while PCM is a related sampling based non-intrusive (and hence 
treats the system ([T]) as a black box) way of implementing gPC. 

2.1.1 Generalized Polynomial Chaos 

In the finite dimensional random space Ps defined in the gPC expansion 
seeks to approximate a random process via orthogonal polynomials of ran- 
dom variables. Let us define one-dimensional orthogonal polynomial space 
associated with each random variable ^/t , /c = 1 , • • • , p as 

VT^^A = {^.rj^^R:ve span{Vi(Cfc)}fio}' (7) 

where, {V'j(Cfc)}f=o denotes the orthonormal polynomial basis from the so 
called Wiener- Askey polynomial chaos . The Askey scheme of polynomials 
contains various classes of orthogonal polynomials such that their associated 
weighting functions coincide with probability density function of the under- 
lying random variables. The corresponding P-variate orthogonal polynomial 
space in Ps is defined as 

W^'^ = W''"^', (8) 
|d|eP 

where the tensor product is over all possible combinations of the multi-index 
d = (di, d2, • • • , d|s|) G nI^I in set P, 

|s| _ 

P = {d G nI^I : |d| = ^ < P and di < Pi} (9) 

i=l 

and, P = (Pi, • • • , P|s|)"^ ^ ^^^^ is vector of integers which restricts the max- 
imum order of expansion of i-th variable to be Pi, and P = maxj Pj. Thus, 
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W^'^ is the space of N-variate orthonormal polynomials of total degree at 
most P with additional constraints on individual degrees of polynomials, 
and its basis functions ^'^ ' {^) satisfy 

/ ^f'''{0^f^iC)pJ,{OdC = S^„ (10) 

for all I < i,j < N-£ = dim{W^'^). Note that in standard gPC all ex- 
pansion orders are taken to be identical i.e. Pi,= i-*2 • • • = -P|s| = so 

that dim(H^^'P) = ^"^tl^P' • We will however take advantage of an adaptive 



expansion, a notion which will be fully developed in section 4.3 



The major advantage of applying the gPC is that a random differential 
equation can be transformed into a system of deterministic equations. A 
typical approach is to employ a stochastic Galerkin projection, in which all 
the state variables are expanded in polynomial chaos basis with correspond- 
ing modal coefficients {aik(t)), as 

Xk{t,0^xf''{t,0 = ^a^kmf'''{0, k = l,---,n, (11) 

i=l 

where, the sum has been truncated to a finite order. Substituting, these 
expansions in Eq. ([T]), and using the orthogonality property of polynomial 
chaos ( 10 ) , we obtain for A: = 1 , • • • , n, j = 1 , • • • , Nj] , 

ajk = Fjk{^,t), (12) 

a system of deterministic ODEs describing the evolution of the modal coef- 
ficients, with initial conditions 

ajk{o)= [ ^fc(o,e)^f''(Ops(e)de, (13) 

and a = (an, • • • , a^j^i, • • • , ai„, • • • , a^^n)'^, 

F,fe(a,t)= / /fc(x^'^(e,t),6,t)^f''(6PE(e)d^, (14) 

with x^'^(t, ,^) = {xf'^{t,^), • • • , Xn'^{t, ^)). This system can be solved with 
any numerical method dealing with initial- value problems, e.g., the Runge- 
Kutta method. Similarly, the observable can be expanded in gPC basis, 
as 

i=l 
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where, 

bMt) = I ^'^mT'^iOpmdc (16) 



with k = 1, ■ ■ ■ ,d. Hence, once the solution to the system (12) has been 
obtained, the coefficients bjk can be approximated as 

bjk- [ 5fc(x^'^(t,0)^f''(0Ps(e)de. (17) 

Such Galerkin procedures have been used extensively in the literature. In 
many instances Galerkin projection may not be possible due to unavailability 
of direct access to the system equations ([T]). In many other instances such 
intrusive methods are not feasible even in cases when the system equations 
are available, because of the cost of deriving and implementing a Galerkin 
system within available computational tools. To circumvent this difficulty, 
probabilistic collocation method has been developed. 

2.1.2 Probabilistic Collocation Method 

PCM is a non-intrusive approach to solving stochastic random processes 
with the gPC. Instead of projecting each state variable onto the polynomial 



chaos basis, the collocation approach evaluates the integrals of form (16) by 
evaluating integrand at the roots of the appropriate basis polynomials j5]. 
Given a ID probability density function Pj{(,j), the PCM based on Gauss 
quadrature rule, approximates an integral of a function g with respect to 
density PjHj)-, as follows 

1 "^'j 

/ 9{ij)p{ij)dij ^l^ij[9\ = '^wi^kg{n^k), j = i, (18) 
•^"^ k=l 

where, r/ ^ G C;. is the set of Gauss collocation points with associated 
weights wi-k, Ij is the accuracy level of quadrature formula, and mi- is the 
number of quadrature points corresponding to this accuracy level. Building 
on the ID quadrature formula, the full grid PCM leads to following cubature 
rule, 

-1 ri /■! 



Ah,-- - Jp,p)[g] = {Ui,(^Ui,---Ui^)[g] 

niip 

j;--- j;u;ij<7(rij), (19) 
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where, ry = {ri^j^ , • • • , n^j^ ) , with 1 = (/i , • • • , /p) , and j = (ji , • • • , jp) and 
wij = wi-^j^ ■ ■ ■ wi^j^ . To compute p) we need to evaluate the function on 
the full collocation grid C(l,p) which is given by tensor product of ID grids 

C{l,p) = Ci,x---xCi^, (20) 

with a total number of collocation points being Q = 11^=1 "^'j • ^^^^ 
framework, therefore, for any t, the approximations to the model coefficients 



ajk{t) (see Eq. (11)) and bjk{t) (see Eq. (15)) can be obtained as 



ajkit) = [ x^'^(t,e)*f''(Ops(Orfe 



^ • • • ^ -"^ij^f ■^(rij)xfc(t,rij), 

ii=i ip=i 



(21) 



with similar expression for hjk{t). Note to compute summations arising in 
(21), the solution x(t,rij) of the system ([T]) is required for each collocation 
point rij in the full collocation grid C{\,p). Thus, simplicity of colloca- 
tion framework only requires repeated runs of deterministic solvers, without 
explicitly requiring the projection step in gPC. 

If we choose the same order of collocation points in each dimension, i.e. 
mi^ = rrii^, ■ ■ ■ = mi^ = I, the total number of points is Q = P, and the 
computational cost increases rather steeply with the number of uncertain 
parameters p. Hence, for large systems (n ^ 1) with large number of 
uncertain parameters (p 3> 1), PCM becomes computationally intensive. 
As discussed in the introduction, alleviating this curse of dimensionality 
is an active area of current research p]. In this paper we propose a new 
uncertainty quantification approach which exploits the underlying network 
structure and dynamics to overcome the dimensionality curse associated 
with PCM. The key methodologies for accomplishing this are the graph 
decomposition and waveform relaxation, which are discussed in subsequent 
sections. 



3 Graph Decomposition and Waveform Relaxation 
3.1 Waveform Relaxation 

In this section we describe the basic mathematical concept of the Waveform 
Relaxation (WR) method for iteratively solving the system of differential 
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equations of the form ([T]). For purposes of discussion here, we fix the pa- 
rameter values ^ in the system ([T]) to fixed mean values. The general struc- 
ture of a WR algorithm for analyzing system ([T]) over a given time interval 
[0, T] consists of two major steps: assignment partitioning process and the 
relaxation process [IHl [l5] . 

Assignment-partitioning process: Let M = {1, • • • , n} be the set of state 
indices, and Ci,i = 1, ■ ■ ■ , m be a partition of M such that 

m 

M=\JCi, Cif|C, =0,Vi/j. (22) 

i=l 

We shall represent hy D : A/"— )• Ai = {1,2,-- - ,m} a map which assigns 
the state index to its partition index, i.e. = j where, j is such that 

i £ Cj. Without loss of generality, we can rewrite Eq. [Tjafter the assignment- 
partitioning process as: 



(23) 

(24) 
(25) 

(26) 
(27) 

are the subvectors assigned to the i—th partitioned subsystem, such that 
jieCi, k = 1,- ■■ ,Mi = \Ci\ and 

d.(t)^iyl,...,ylf, (28) 

is a decoupling vector, with jk G Mi and k = I, - ■ ■ ,Ni = \Mi\. Here, Mi 
is the set of indices of the partitions (or subsystems) with which the i—th 
partition (or subsystem) interacts, and is given by 

Ar, = M\^{V{Ar.n), (29) 





yi = 


Fi(yi,di(t),Ai,t) 




Ym = 


f ?n(yr?i) 


djn{t),km,t) 


where, for each i = 1, 


■■■ ,m, 










- (/;!'••■ 






yi = 






with initial condition 










yoi = 






and 
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where, = {j £ M : |^ = 0} and denotes the image of the map V. 

Relaxation Process: The relaxation process is an iterative procedure, 
with fohowing steps 

• Step 1: (InitiaHzation of the relaxation process) Set 1 = 1 and guess 
an initial waveform {y['(t) : t G [0 T]} such that y^(0) = yoi, i = 
1 • • • , m. 

• Step 2: (Analyzing the decomposed system at the I-th WR iteration) 
For each i = 1, • • • , m, set 

df(t)=-((yj-^)^,---,(yl-Tf, (30) 

and solve the subsystem 

yf = F,(yf,df(t),A„t), (31) 

over the interval [0,rs] with initial condition yf (0) = yoj, to obtain 
{y'{t):te[0,Ts]}. 

• Step 3 Set 1 = 1+1 and go to step 2 until satisfactory convergence is 
achieved.. 

The general conditions for convergence of WR for a system of differential 
algebraic equations (DAEs) can be found in |18tll7j. Here, we recall a result 
from [17j specializing it for a system of differential equations. 

Proposition 3.1. Convergence of WR for ODE's (see flS^ for proof): 
Given that the system |ip is Lipschitz ( condition then for any initial 
piecewise continuous waveform {yi'(t) : t G [0,Ts]} such that y-'(O) = yoi 
(see definition {2(^), i = 1- ■ ■ ,m, the WR algorithm converges to the solu- 
tion of ([Ip with initial condition xq. 

A more intuitive analysis of error at each waveform iteration is described 
in [T7]. Let y be the exact solution of the differential equation (23 ) and define 
Ej to be the error of the I-th. iterate, that is 

Ej{t)=y'{t)-y{t). (32) 

As shown in |17j . the error \Ej\ on the interval [0, T] is bounded as follows 

\Ei{t)\ < ^^^\Eo{t)\, (33) 
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with C = e^'^ . Here C and 77 are related to the Lipschitz constants of the 
waveform relaxation operator |T7]. It is important to note here that the 11 in 
the denominator dominates the error. Thus, with enough iterations one can 



make the error fall below any desired threshold. It is also evident from Eq. 33 
that the error of standard waveform relaxation crucially depends on T. The 
longer the time interval, the greater is the number of iterations needed to 
bound the error below a desired tolerance. Based on this observation, a 
novel "adaptive" version of waveform relaxation has been developed in 117] , 
which we refer to as adaptive waveform relaxation (AWR). The idea here 
is to perform waveform relaxation in "windows" of time that are picked so 



as to reduce / in Eq. 33, Specifically, one can pick small time intervals for 
computation when the solution to ([T]) changes significantly (implying EQ{t) is 
large) and pick large intervals when the solution changes little (consequently 
£'o(t) is small). The solution from one time interval is extrapolated to the 
next using a standard extrapolation formula [19] and the initial error is 
estimated using, 

^/+i,o(sJ (T+I)! ^ ^ ' 

where, 

u{t) = {t-tQ){t-ti)...{t-ti). (35) 

Here to, ti, • . • , are points through which one passes the extrapolating poly- 
nomial [IT]. Note that, (f)^^^ = ^ is the l-ih. derivative of the waveform 
relaxation operator (/> with respect to t (see. [T^ [T7]). 

The algorithm to compute the length of the time windows is as follows: 
Adaptive Waveform Relaxation: To compute the time interval for AT/+1, 
execute the following steps: 

1. Set ATz+i = 2Arf and 5 = ^ATi. 



2. Evaluate EiJ^ifi{Tj + ATf+i) using Eqn. 34 to estimate Hi^j+i^oll and 



compute with the aid of the following equation, 

\\Ei-,,,r\\ = |^^^ll^m,o||. (36) 

3. If > £ and AT/+1 > ^ATj, set ATf+i = AT/+1 - 6 and 

repeat step 2. 
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We define the minimal window length to be AT/4_i = ^T. The above proce- 
dure gives a sequence of time intervals [0,Ti], [Ti,T2], . . . , [T^-i,Tiy], where 
Tu = T, on which WR is performed (as described in relaxation process) with 
an initial "guess" waveform provided by an extrapolation of the solution on 
the previous interval [IT]. AWR is found to accelerate simulations by or- 
ders of magnitude over traditional WR |17] . In this work, we propose to use 
AWR for simulating the set of differential equations that arise from intrusive 
polynomial chaos. As mentioned before, the curse of dimensionality gives 
rise to a combinatorial number of equations [1] making AWR particularly 
attractive. 

While the convergence of WR or AWR is guaranteed irrespective of how 
the system is decomposed in the assignment-partitioning step, the rate of 
convergence depends on the decomposition |17j . For a given nonlinear sys- 
tem, determining a decomposition that leads to an optimal rate of AWR 
convergence is an NP-complete problem [T7j. Ideally, to minimize the num- 
ber of iterations required for convergence, one would like to place strongly 
interacting equations/ variables on a single processor, with weak interactions 
between the variables or equations on different processors. We show in |17j . 
that spectral clustering [T3] along with horizontal vertical decomposition [T2] 
is a good heuristic for decomposing systems for fast convergence in WR and 
AWR. For this task, we now discuss a novel decentralized spectral clustering 
approach [16J that when coupled with AWR [17j provides a powerful tool 
for simulating large dynamical systems. 

3.2 Graph Decomposition 

The problem of partitioning the system of equations ([T| into subsystems 
based on how they interact or are coupled to each other, can be formulated 
as a graph decomposition problem. Given the set of states xi, • • • and 
some notion of dependence Wij > 0, i = 1, • • • ,n, j = 1, ■ ■ ■ , n between pairs 
of states, an undirected graph G = {V, E) can be constructed. The vertex 
set y = {1, . . . , n} in this graph represent the states Xj, and the edge set is 
E (^V X V ^ where a weight Wij > is associated with each edge (i, j) G -E, 
and W = [wij\ is the n x n weighted adjacency matrix of Q. In order to 
quantify coupling strength idij between nodes or states, we propose to use 

Wij = ^[\Jij\ + \Jji\], (37) 
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where, J = [^^ //q""^"^" "^iji^i^'i ^m) , ^m,t)dt], is time average of the Jacobian, 



J(x,e,t) 



a/,(x(t;0,C,i) 



dxj 



(38) 



computed along the solution x(t; ,^) of the system ([T]) for nominal values 
of parameters ^m- Use of system Jacobian for horizontal vertical graph 
decomposition can also be found in |12j . 

We will now discuss, spectral clustering (see [13j a popular graph de- 
composition/clustering approach that allows one to partition a undirected 
graph given its adjacency matrix W . In this method first a (normalized) 
graph Laplacian is constructed as follows [U [201 E] , 

(1 ifi = j 

Lij = < if (^,i) e E (39) 

otherwise , 

or equivalently as, L = I — D~^W where D is the diagonal matrix with the 
row sums of W. The clustering assignment/decomposition is then obtained 
by computing the eigenvectors/eigenvalues of L. In particular, one uses the 
signs of the components of the second and higher eigenvectors to partition 
the nodes in the graph into clusters |I3] . Traditionally, one can use standard 
matrix algorithms for eigenvectors/eigenvalues computation |22| . However, 
as the size of the dynamical system or network (and thus corresponding 
adjacency matrix) increases, the execution of these standard algorithms be- 
comes infeasible on monolithic computing devices. To address this issue, 
distributed eigenvalue/eigenvector computation methods have been devel- 
oped, see for example [23] . 

In [16j, a wave equation based distributed algorithm to partition large 
graphs has been developed which computes the partitions without construct- 
ing the entire adjacency matrix W of the graph |13j . In this method one 
"hears" clusters in the graph by computing the frequencies (using Fast 
Fourier Transform (FFT) locally at each node) at which the graph "res- 
onates". In particular, one can show that these "resonant frequencies" are 



related to the eigenvalues of the graph Laplacian L ( 39 ) and the coefficients 
of FFT expansion are the components of the eigenvectors. Infact, the al- 
gorithm is provably equivalent to the standard spectral clustering [13], for 
details see [16] . 

The steps of the wave equation based clustering algorithm are as follows. 
One starts by writing the local update equation at node i based on the 
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discretized wave equation, 



Ui{t) = 2ui{t - 1) - Ui{t - 2) - ^ LijUj{t - 1) (40) 

where Ui{t) is the value of u at node i at time t and Ljj are the local 
entries of the graph Laplacian. At each node i, Uj(0) = 1) is set to a 
random number on the interval [0, 1]. One then updates the value of Ui using 



Eqn. (40) until t = T^ax (for a discussion on how to pick T^ax see [16]). 
Note that, Ui{t) is a scalar quantity and one only needs nearest neighbor 
information in Eqn. ( |40[ ) to compute it. One then performs a local FFT on 

[ui{l), ,Ui{Tmax)] and then assigns the coefficients of the peaks of the 

FFT to Vij. Here the frequency of the j-th peak is related to Xj, the j-th 
eigenvalue of the graph Laplacian L, and Vij is the i-th component of the 
j-th eigenvector. 

In |16j . it has been shown that for the wave speed c < \/2 in Eqn. 



(40) above wave equation based iterative algorithm is stable and converges. 
Moreover, the algorithm converges in 0{^/T) steps, where r is the mixing 
time of the Markov Chain [16j associated with the graph G. The competing 
state-of-the-art algorithm [23] converges in 0(r(log(n)^). For large graphs 
or datasets, 0{^/T) convergence is shown to provide orders of magnitude 
improvement over algorithms that converge in 0(r(log(n)^). For detailed 
analysis and derivations related to the algorithm, we refer the reader to |16| . 



4 Scalable Uncertainty Quantification Approach 

In this section we discuss how gPC and PCM can be integrated with WR 
scheme extending it to a probabilistic setting. As mentioned earlier we re- 
fer to this iterative UQ approach as PWR. Figure [T] shows the schematic 
of PWR framework. In the intrusive PWR, the subsystems obtained from 
decomposing the original system are used to impose a decomposition on sys- 
tem obtained by Galerkin projection based on the gPC expansion. Further 
the weak interactions are used to discard terms which are expected to be 
insignificant in the gPC expansion, leading to what we call an Approximate 
Galerkin Projected (AGP) system. We then propose to apply standard or 
adaptive WR on the decomposed AGP system to accelerate the UQ com- 
putation. In the non-intrusive form of PWR, rather than deriving the AGP 
system, one works directly with subsystems obtained from decomposing the 
original system. At each waveform relaxation iteration we propose to apply 
PCM at subsystem level, and use gPC to propagate the uncertainty among 
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Figure 1: Schematic of intrusive (left) and non-intrusive (right) PWR. 



the subsystems. Note that unhke intrusive PWR (where deterministic de- 
couphng vectors or deterministic waveforms are exchanged) , in non-intrusive 
PWR, stochastic decouphng vector or probabihstic waveforms represented 
in gPC basis are exchanged between subsystems at each iteration. 

We first describe the key technical ideas behind intrusive and non-intrusive 



PWR though an illustration on a simple example in section 4.1 These no 



tions are fully generalized later in sections 4.2 4.4 We also prove the con 



vergence of PWR approach (in section 4.5), and in section 4.6 discuss the 



computational gain it offers over standard application of gPC and PCM. 
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4.1 Main Ideas of PWR 



We illustrate the proposed PWR framework through an example of para- 
metric uncertainty in a simple system ([T]). Consider the following coupled 
oscillator system: 

xi = fi{xi,X2,uJi,t) = uji + Ki2sm{xi- X2), 

X2 = f2{xi,X2,X3,UJ2,t) 

= clI2 + -fC2i sin(xi - X2) + -^^23 sin(x3 - X2), (41) 

X3 = f3{x2,X3,UJ3,t) = UJi + Ks2Sm{x2 - X3) 

Here, uji is the uncertain angular frequency of the i*'^ (i = 1,2,3) oscilla- 
tor. We assume that parameters are mutually independent, each having 
probability density pi = piiuji). The coupling matrix K is 

/ K12 \ 
K= \ K21 K23 (42) 
V K32 / 

is assumed deterministic with Kij = 0(e), e ^ 1, so that the three oscillators 
in ([1]) weakly interact with each other, i.e. the subsystem 2 weakly affects 
subsystem 1 and 3, and vice versa. 



4.1.1 Approximate Galerkin projection for the simple example 

In standard gPC, states Xi are expanded in a polynomial chaos basis as 

'(i,a;) = ^a,-,(f)^f^(u;), i = 1,2,3 (43) 



where, ^ ■ ' G W^'^ , the P variate polynomial chaos space formed over 
S = {uJi,uj2,uJ3} and P = {Pi, P2, P3) determines the expansion order(see 



section 2.1.1 for details). Note that in this expansion (43), the system states 
are expanded in terms of all the random variables uj affecting the entire 
system. From the structure of system ([l]) it is clear that the 1** subsystem 
is directly affected by the parameter wi and indirectly by parameter uj2 
through the the state X2- We neglect second order effect of cjs on xi. A 
similar statement holds true for subsystem 3, while subsystem 2 will be 
weakly influenced by ui and UJ3 through states xi and X3, respectively. This 
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structure can be used to simplify the Galerkin projection as follows. For xi 
we consider the gPC expansion over Si = Ai |J A^, 

xf''^'{t,ui,u;2) = X^OjiW^f '^'(^^i,c^2), (44) 

where, 

Ai = {wi}, A^ = {0.2} (45) 

and Pi = {Pii,Pi2)- Note that since uj2 weakly affects xi, the order of 
expansion P12 can be chosen to be smaller compared to Pn. Similarly, 
one can consider simplification for X3^''^^(t,a;i, ws). For X2 following similar 
steps, we define 

A2 = {a;2}, A^ = {wi,W3} (46) 

and P2 = {P21, P22, P23)- By similar argument, one will choose P2i,P23 
much smaller than ^22- We also introduce the following two projections 
associated with the state X2- 

p2.(^S„P.) ^ J2 ^^E„P.^ ^E.P.^ ^E„P,_ (47) 

i=i 

where i = 1,3 and (•,•) is the appropriate inner product on W^'^ (see 
section 4.3 for details). With these expansions, and using standard Galerkin 



projection we obtain the following system of deterministic equations 

5 = F(a,t), (48) 
with appropriate initial conditions, where 

i^ji(a) = / Mxf^^''\r^'Hx2"'^'),^i,t)^f'''''{^)pi^)doo,, 

Fj2{si) = [ f2ix^"''\x^''''\x^"''\L02,mf''''{u;)p{u^)du;, 

F,3(a) = / /3(P=^'3(x^^'^^),xf^'^\a;3,t)*f'''^(^)Ku;)dw, 

and a = (ai, £12,^3)^, with a^ = (aji, • • • ..aiN^^) and F = (Fi,F2,F3)^ with 
Fj = (Fji,--- ,Fijv^ ). We will refer to (48) as an approximate Galerkin 
projected (AGP) system. The notion of aGP in more general setting is 
described in section |4^ 
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4.1.2 Intrusive PWR illustrated on the simple example 



In intrusive PWR, after performing the AGP explicitly, the system (48) is 
decomposed as 

Qij = Fji{sLi,di,t), (49) 

where di = 7^^'^ (212), d2 = (ai,a3) and da = ^^'^(£12) are the decoupling 
vectors (here we overloaded notation for 'P^'*(a2) to imply the coefficients 



in expansion (|47[) ) . Note that the decomposition of system ( 48 ) is based on 



the decomposition of the original system (41). Adaptive or standard WR 



described in section 3.1 , can then be applied to solve the decomposed system 



(49), iteratively. Since, the the system (49) is deterministic, deterministic 
waveforms or deterministic decoupling vectors dj,i = 1,2,3 are exchanged 
in each WR iteration (see Figure [T] for an illustration) . 

4.1.3 Non-Intrusive PWR illustrated on the simple example 

In the non- intrusive form of PWR, rather than deriving the AGP, one works 
directly with subsystems obtained from decomposing the original system. 
The main idea here is to apply PCM at subsystem level at each PWR iter- 
ation, use gPC to represent the probabilistic waveforms and iterate among 
subsystems using these waveforms. Recall that in standard PCM approach 



(2.1.2), the coefficients a^(t) are obtained by calculating integral 



The integral above is typically calculated by using a quadrature formula and 
repeatedly solving the i*^ subsystem over an appropriate collocation grid 
C*(Sj) = C*(Aj) xC*(A?), where, C*(Ai) is the collocation grid corresponding 
to parameters Aj (and let Ig be the number of grid points for each random 
parameter in Aj), C*(A?) is the collocation grid corresponding to parameters 
A^ (and let Ic be the number of grid points for each random parameter in A? 
). Since, the behavior of i^^ subsystem is weakly affected by the parameters 
A?, we can take a sparser grid in A| dimension, i.e. Ic < h, as we took 



lower order expansion for these random variables in section 4.1.1 Below we 
outline key steps in non intrusive PWR: 

• Step 1: (Initialization of the relaxation process with no coupling ef- 
fect): Set 1 = 1, guess an initial waveform x^{t) consistent with initial 
condition. Set d\ = X2, = {xi,x'^), dg = X2, and solve 



X. 



l = fi{x},d}{t),ooi,t), (51) 



18 



with an initial condition xj{0) = x^{0) on a collocation grid C*(Aj). 
Determine the gPC expansion xf"^"^(t, •) by computing the expansion 



coefficients from the quadrature formula (50) 



Step 2: (Initialization of the relaxation process, incorporating first 
level of coupling effect): Set 1 = 2 and let = x^^'"^^'^, d2 = 
{x^^'^^'^ , x^^'^^'^), d| = X2^'^^'^ be the stochastic decoupling vectors. 
Solve 

x^ = fi{xld^,{t,-),uji,t), (52) 

over a collocation grid C*($]j) to obtain x- " (t, •). From now on 
we shall denote the solution of the i*'^ subsystem at iteration by 

Step 3 (Analyzing the decomposed system at the I-th iteration^: Set 
the decoupling vectors, d( = V'^'^{x^^'^^'^~^), d2 = Xg ^'■^^'^~^) 

= p2,3(^S2,P2,/-l^ g^j^g 

xi = fi{xi,di{t,-),u:i,t), (53) 

over a collocation grid C*(Sj) and obtain the expansion x^ " " (t, •). 

Step 4 (Iteration) Set / = / + 1 and go to step 5 until satisfactory 
convergence has been achieved. 



Note that in above non-intrusive PWR, the decoupling vectors are stochastic 
and so at each iteration probabilistic waveforms are exchanged between sub- 
systems ((see Figure [T] for an illustration). We next generalize the intrusive 
and non-intrusive PWR introduced above, in the forthcoming sections. 

4.2 Decomposition of Galerkin Projected System 

We begin by revisiting the complete Galerkin system ( [T2| ). To apply WR, 



recall that the first step is the assignment-partitioning (see section 3.1) 
There are two possible approaches for partitioning the complete Galerkin 
system. One can first split the original dynamical system ([T]), and then use 



this decomposition to partition the complete Galerkin projection (12) by 



assigning the model coefficients in (11 ) for each state to the cluster to which 



state is assigned to while decomposing system ([T]). As previously explained 
in section [3^ the partitioning is performed by representing the dynamical 
system ([T]) as a graph with the symmetrized time averaged Jacobian (37) 



as the weighted adjacency matrix. One can then apply the wave equation 



based decentralized clustering algorithm outlined in section 3.2 
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Alternatively, one can perform this decomposition directly on the com- 



for the resulting system ( 12 ) be 



plete Galerkin projection (12). Let the symmetrized time averaged Jacobian 



where, J 



[jr J^°~^'^" Jij{a{t),t)dt], is time average of the Jacobian, 



J (a, t) 



5Fifc(a(t)) 



da 



ik 



computed along the solution a(t) of the system (12). This gives, 

r 9A(x^(e,t),6,t), 



J (a,t) 



dajk 



Taylor expanding fk{x{^,t),(,k,t) locally, gives, 

j'(a,t) ^ J{a,t). 



(54) 



(55) 



(56) 



(57) 



Thus, one expects to get similar results by performing clustering on the 
original system (in ([T])) to that obtained based on complete Galerkin sys- 
tem (12). Since the dimensionality of system ([T]) is much lower than that 
of system (12), the first decomposition is less computationally challenging 



than the latter. In this work, we will use the original system to determine 
the decomposition and use that to impose the partition of the Galerkin pro- 



jection. Given the decomposition of system (12), one can use the WR or its 



adaptive form to simulate the system in a parallel fashion. However, before 
doing this one can further exploit the weak interaction between subsystems 
to reduce the dimensionality of the complete Galerkin system, as described 
in section 4.1.1[ We next describe this approximate Galerkin projection in 
a more general setting. 



4.3 



Approximate Galerkin Projection and Intrusive Proba- 
bilistic Waveform Relaxation 



Recall that in the gPC expansion (11), all the system states are expanded 
in terms of random variables affecting the entire system. However, the 



i—th. subsystem in the decomposition (see Eq. (23)) is directly affected by 
the parameters A, (see definition (27)) and indirectly by other parameters 
through the decoupling vector (see definition ([28]) ) . We shall denote by 



A? 



■J' 



(58) 
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the set of parameters which indirectly affect the i— th subsystem through 
the immediate neighbor interactions and by 

Si = A, U A^, (59) 

the set of parameters that directly and indirectly (through nearest neigh- 
bor interaction) affect the i— th subsystem. Under the hypothesis that the 
i—th subsystem is dynamically weakly coupled with its nearest neighbors, 
uncertainty in parameters A* will weakly influence the states in i—th sub- 
system through the decoupling vector, while the uncertainty in parameters 
S \ Sj can be neglected. To capture this effect, consider a P-variate space 
(analogous to the P-variate space introduced in the section 2.1.1) 

W^'^ = (g) W''^''^'^^, (60) 

|d|GP 

formed over any random parameter subset A = {Cji)Ci2) ' ' ' ^Cin} ^- We 
shall denote the basis elements of W^'^ by = I,-- - ,N/^, where 

TVa = dim{W^'^). Note that for any Ai C A2 C S, 

where, W'^ = {0} is the P— variate space corresponding to the empty set. 
Also, recall that Pi, P2 are vectors which control the expansion order in gPC 
expansion. The inner product on W^'^ induces an inner product on W^'^ 
as follows 

<Xi(0,X2(0 >A= / Xi{0X2{C)PA{0dt (62) 

for any Xi(S^), X2{^) G W^'^ . Using this inner product, we introduce a 
projection operator 

p^A2 . ^A2,P2 ^ ^Ai,Pi^ (63) 

such that for any X{C) G W^^'^^ 

PrtSli^m = E < >A ^t^'^'iO- (64) 

i=l 

With respect to the given decomposition D imposed on the system (see 
section 3.1), we define a projection operator V^'^ indexed by subsystem i 



and state xj 



V 



Ao(,)U{Aj,(^)nA.), 



Prf^ if P(j)/i,MnAA^( 
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Remark 4.1. For any subsystem i, since the parameters A? weakly affect it, 
we can adaptively select component of vector Pi = (Pn,--- ,-PjjSi|) -^o that 
a lower order expansion is used in components corresponding to random 
variables in A?. 

For any subsystem i and a vector x^'^(^, t) = t), - ■ ■ , t)) 

(where ' {C,t) is standard gPC expansion (11)), we associate a vector 
valued projection operator as follows 

r(x^'^) = (r-'^Hx^;^),... (65) 

In terms of these operators, for any state Xk, an approximate Galerkin pro- 
jected equation is defined as, 



|r''=[xf'^(e,t)] = /fc(p*(x^'^(e,t)),6,t), (66) 



where, i = T){k) is the index of the subsystem to which the state k belongs. 
More, precisely the above system can be expressed as: 

a]k = F]k{a-,t), (67) 



where, a = (af/^^ , • • • , a^[,^^ , • • • , afj"^ , • • • , a^^f' 



F)k = [ /fc(^'(x^'''(e,i)),Cfc,i)ff''''(e)ps.(0^e, 
and, j = 1, - ■ ■ , Nj]^, k = 1, - ■ ■ , n with i = V{k). Let 



then the system (67) can be compactly written as 

n = F(a,t), (69) 



with appropriate initial condition (see expression 13) and will be referred 
to as the approximate Galerkin projected (AGP). Using this generalization 
of AGP system, it is straightforward to generalize the intrusive PWR intro- 
duced in section 14.1.21 
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4.3.1 Intrusive PWR Algorithm 



In the intrusive PWR, one applies tlie WR to AGP system. As per discussion 
in section 4.2 the decomposition D on the original system ([T]) is used to 
imposes a decomposition on the system (69) leading to 



for i = 1, 



, m, where 



ai = Fj(ai,di,t), 



(70) 



(71) 



(72) 



ki S Cj, and dj = (a^, • • • ,aj^ is the decoupling vector (recall notation 
from section [3^1] ) . One then follows the procedure for waveform relaxation 
or its adaptive version, as described in section [3?T| Adaptive WR can lead 
to a significant increase in convergence of WR as demonstrated in |17| , and 
would be illustrated later in the section [H 



As discussed in section 2.1.1 , the projection step (66) can be very tedious 



and in some cases not possible. Hence, applying waveform relaxation directly 



to the system (70) may not be practical in many instances. In the next 



section, we describe an alternative non-intrusive approach using probabilistic 



collocation, which does not require the projection step (66) explicitly. 



4.4 Non-Intrusive Probabilistic Waveform Relaxation 



In terms of the projection operator (65), we can rewrite each subsystem in 



(23) as 



Yi = Fi{yi,V\di{t,-)),Ai,t), (73) 
where, dj(t, •) is the stochastic decoupling vector or probabilistic waveform, 

d.(V) = ((y^'^^)^,---,(yS"'''^"OT. (74) 



where, we have explicitly indicated the dependence on the parameters ( see 



definition (28)). Here for any i = 1, • • • ,m, y 
with 



Jl J Mi 



x^'^^'iCt) = (t)^'^'^(6 = r^'^k^:''), 



m=l 



Ik 



(75) 
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By definition the coefficients a^^i (i) in above expansion satisfy the system 
67|). These coefficients can be obtained by using the quadrature formula 



21| , by repeatedly solving the system ( 73 ) over an appropriate collocation 
grid C(l,nj) 

C(l, m + nl) = C(o, m) X C(m, <), (76) 

where, 1 = (o,m), C(o, n^) = Cl x ■ ■ ■ x Cl , is the collocation grid 
corresponding to parameters Aj, with rii = |Aj|, o = (oi, • • • , OnJ, and 
C(m, n?) = X • • • X ^, is the collocation grid corresponding to pa- 

rameters A?, with n? = | A?| and m = (mi, • • • , nin'^). For simplicity we take 
oi = ■ ■ ■ = Om = Is and mi = • • • = mn^^ = /c for i = 1 • • • , m. Since, the 
behavior of i— th subsystem is weakly affected by the parameters A? through 



the decoupling vector, then consistent with remark (4.1) we can take 



Ic < Is, (77) 

leading to an adaptive collocation grid for each subsystem. With this, we are 
ready to generalize the non-intrusive PWR approach introduced in section 
11X31 

4.4.1 Non-Intrusive PWR Algorithm 

• Step 1: Apply graph decomposition (see section [3] for details) to iden- 
tify weakly interacting subsystems in the system ([T]). 

• Step 2 (Assignment-partitioning process): Partition ([T| into m sub- 
systems (obtained in Step I) leading to system of equations given by 



(23). Obtain, Aj, A? and Sj for each subsystem, i = 1, • • • , m. Choose 



the parameters I si, Id, Pi- 

• Step 3: (Initialization of the relaxation process with no coupling ef- 
fect): Set 1 = 1 and , guess an initial waveform {y^(t) : t G [0,Ts]} 
for each i = 1, • • • , m consistent with initial condition (see Step 1 in 



relaxation process described in section 3.1). Set 

dj{t) = {yj,{t),--- ,yj^S*)), (78) 
i = 1, • • • , m, and solve for {yf "^"""^(t), t G [0, Ts]} using 

y} = F\y},d}{t),A„t), (79) 

with an initial condition yj^(O) = y^(0) on a collocation grid C(o, n,). 
Determine the gPC expansion yf''''^"^{t, •) over P— variate polynomial 
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space W^^'^"^ by computing the expansion coefficients from tlie quadra- 



ture formula (21). 



• Step 4: (Initialization of the relaxation process, incorporating first 
level of coupling effect): Set 1 = 2 and for each i = 1, • • • , m, set 

d?(V) = (y^'''-\t, •),•••, y^;:''''^"-^^)), (80) 

and solve for {yf{t),t G [0,Ts]} from 

y2=F(y?,df(t,.),A.,t), (81) 

with an initial condition yf{0) = YiiO), over a collocation grid C{1, ni+ 
nf). Obtain the expansion " " (t, •) using ( |75[ ). From now on we 
shall denote the solution vector of the i—th subsystem at /— th itera- 
tion by yf 

• Step 5 (Analyzing the decomposed system at the I-th iteration): For 
each i = 1, • • • , m, set 

<l.' = {yJ'''-"-",....yJ-'''"''"-"), (82) 

and solve for {y^{t) : t G [0,Ts]} from 

yf = r(yf,P*(df(t,.)),A„t), (83) 
with initial condition y^ (0) = yf(0) over a collocation grid C(l, nj+n?). 



Obtain the expansion yf"^"\t, •) using the expansions (75). 



• Step 6 (Iteration) Set 1 = 1 + 1 and go to step 5 until satisfactory 
convergence has been achieved. 

4.5 Convergence of PWR 

Below we prove that the iterative PWR approach converges. The proof is 
based on showing that the AGP system is Lipschitz if the orginal systems 
is Lipschitz (see condition [2]) , and then invoking standard WR convergence 



result (3.1 ). 



Proposition 4.2. Convergence of PWR: The intrusive and non-intrusive 



PWR algorithms described in sections 4-3- 1 and J^.^-.l, respectively converge. 
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Proof: We prove the result for intrusive PWR. By construction, since 
non- intrusive PWR algorithm solves the AGP system ( 69 ) in a different way, 



the convergence result holds for non-intrusive PWR as well. Consider the 

-IV{1) _lV(n) 



AGP system (legl) and let 



ai = (a 



ic(i) 

11 5 ■ ■ ■ ; ' 



and 



_2D(1) 



Let for a given k = 1, ■ ■ ■ ,n, i = T){k), then 



2V{n) 



-2V{n) 



(84) 
(85) 

(86) 



and P^(x''^'^) = (P^'i(x'{^'^),--- ,P^'"(xk^'^)), for I = 1,2 and for sim- 
plifying notation we have dropped subscripts on P vectors. Then for each 
k = !,■■■ ,n,i = V{k), j = !,■■■ ,iVs^, 

||F}fc(a2) -F}fc(ai)|| 

(/fc(r (x^'^'2), e, t) - A(pxx^'^'i), e, t)) 

j 



/IV 1^ \ III, I 

m=l p=l 

x|^f'''(0ks(0rfe 



m=l p=l 



*2'(m)|C l.fCm) _ 2,D(m)^ 



(87) 



where, 



For a given i = 1, - ■ ■ ,n and j = 1 • • • , Nj^^. , let 



'X>(n) 



«7 L jl iNt,^,,, .71 



(1) 
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and define 



then, 

























Lfi. 



|F(a2)-F(ai)|| < L||a2-ai| 



(90) 



where, L = \ \L^\ \ is a matrix norm of . Hence, the system ( 69 ) is Liptchitz. 
Thus, given the original system ([T]) is Lipschitz (condition Q), the approx- 
imate system ( |69[ ) is also Lipschitz as shown above ( |90[ ). Hence, by propo- 
sition 3.1 (adaptation of Theorem 5.3 in ([E])), we conclude that PWR 
converges. 

The question of how the decomposition of the system and the choice of 
the PWR algorithm parameters P, Ig, h influence: 1) the rate of convergence 
of PWR, and 2) the approximation error (due to the truncation introduced 



in the AGP system (69), the use of adaptive collocation grid i.e. condition 
77 and computation of the modal coefficients by the quadrature formula), 
needs to be further investigated. 



4.6 Scalability of PWR 

The scalability of non-intrusive PWR relative to full grid PCM is shown 
in Figure 10, where the ratio IZp/Tti indicates the computation gain over 
standard full grid approach applied to the system ([l]) as a whole. Here 
IZp = is the number of deterministic runs of the complete system ([T]), 
which comprises of m subsystems each with = 1, • • • ,m uncertain pa- 
rameters, such that p = ^™ 1 Pi and / denotes the level of full grid. Similarly, 
TZi = 1+Y1T=1 ^?'+-^max(Z]i^i is the total computational effort 

with PWR algorithm, where /max is the number of PWR iterations. Clearly, 
the advantage of PWR becomes evident as the number of subsystems m and 
parameters in the network increases. Moreover PWR is inherently paralleliz- 
able. 



5 Example Problems 

In this section we illustrate intrusive and non-intrusive PWR through several 
examples of linear and nonlinear networked systems with increasing number 
of uncertain parameters. While most examples are of ODE's, we also give 
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Figure 2: Scalability of PWR algorithm, when implemented with full grid 
collocation as the subsystem level UQ method, with = 5, = 1, • • • , m, 
I = Is = 5, Ic = 3 and /max = 10, 50, 100. The computational gain of PWR 
becomes insensitive to /max, as the number of subsystems m increase. 



an example application of PWR to an algebraic system. This illustrates 
how in principle one can extend application of PWR to DAEs, just like WR 
approach extends to DAEs |15j). Through some examples we study how 
the strength of interaction between subsystems affects the convergence rate 
and the approximation error of PWR. In one of the examples related to 
building model, we also show how time-varying uncertainty can be incorpo- 
rated into standard UQ framework by using Karhunen Loeve expansion. In 
all the examples, we compare solution accuracy of PWR with other UQ ap- 
proaches (e.g. Monte Carlo and Quasi Monte Carlo methods), and wherever 
appropriate mention computation gain offered by PWR over the standard 
apphcation of gPC and PCM. 

5.1 Stability Problem 

We first consider a simple system, with two states (xi,X2) G M^, 

xi = axi+cx2 — vi, (91) 
X2 = cx\ + hx\ — V2-, (92) 

where, c,vi,V2 are fixed parameters, and a,b are uncertain with Gaus- 
sian distribution Q and tolerance 20% (i.e. a = 0.2//, where is the 
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mean and a is the standard deviation of ^). The parameter c determines 
the couphng strength between two subsystems described by the two equa- 
tions. The output of interest is the stabihty of the system, which is deter- 
mined by Am, the maximum eigenvalue of the Jacobian J{xiq, X20; a, ^, c) = 
2aa;io 2cx2o 



2cx 



10 



2bx 



20 



, where, xio,X2o is the equilibrium point satisfying 



QxIq + CX20 - vi 

cx\q + bx^ 



0, 
0. 



(93) 



Figure [s] shows the UQ results obtained by applying PWR (with Ig = 5, 
/c = 3 and Pi = (5,3),P2 = (3,5)) to iteratively solve the algebraic system 



(93). We make comparison with the true (to imply more accurate result 



obtained by solving the complete system (93)) distribution of Am obtained 
by using a full collocation grid on the parameter space (a, b) with la = 5,lb = 
5,P = (5,5). PWR converges to the true mean and variance as shown in 
the left and right panel of the Figure [3] for two different values of c. As the 
coupling strength c increases (see right panel in Figure [s]) , the number of 
iterations required for the convergence increases, as expected. 



c=0.10 



c=2.80 



6.35 





#lter 



Figure 3: Left Panel: Convergence of mean (//) and variance (u) of Am for 
c = 0.1. Right panel: Convergence of mean and variance of Am for c = 2.8. 
Black line indicates the true values. 
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5.2 Building Example 



For energy consumption computation, a building can be represented in terms 
of a reduced order thermal network model of the form [24j , 



where, T E i?" is a vector comprising of internal zone air temperatures, 
and internal and external wall temperatures; ^(u(t);^) is the time depen- 
dent matrix with ^ being parameters, u(i) is control input vector (compris- 
ing of zone supply flow rate and supply temperature), and vectors Qe = 
{Tamb{t),Qs{t))'^ represent the external (outside air temperature and solar 
radiation), and Qi is the internal (occupant) load disturbances. We consider 
the problem of computing uncertainty in building energy consumption due 
to uncertainty in building thermal related properties and uncertain distur- 
bance loads. These uncertainties can be categorized into: (i) static paramet- 
ric uncertainty which include parameters such as wall thermal conductivity 
and thermal capacitance, heat transfer coefficient, window thermal resis- 
tance etc.; and (ii) time varying uncertainties which include the external 
and internal load disturbances. 

Recall, that the traditional UQ approaches and PWR which builds on 
them, can only deal parametric uncertainty. To account for time varying 
uncertain processes, we employ Karhunen Loeve (KL) expansion [25j. The 
KL expansion allows representation of second order stochastic processes as 
a sum of random variables. In this manner, both parametric and time 
varying uncertainties can be treated in terms of random variables. We next 
demonstrate both intrusive and non- intrusive PWR methods. 

5.2.1 Two Zone Example 

We first consider a simplified two zone building model as shown in Fig. |4j 
Here the state T is a 10 dimensional vector comprising of internal wall tem- 
peratures and the internal zone air temperatures, where we have assumed 
that the outer wall surfaces are held at ambient temperature. We also as- 
sume that the ambient temperature and solar load are deterministic fixed 
quantities and there is no internal occupant load. Thus, in computing the 
uncertainty in the zone temperatures, we only consider parametric uncer- 
tainty. Specifically, we assume that the heat transfer coefficient and the 
thermal conductivity of the walls in each zone have standard deviations of 
10% around their nominal values of 3.16W/m?/K and A.GbW / m/ K , respec- 
tively. Thus, locally each zone is affected by two uncertain parameters. 




(94) 
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Figure 4: Diaffram of the two zone tliermal model of a building. To^Ynbij') 
293K in this example. 



with heat transfer coefficient being a common (i.e. same) parameter and 
thermal conductivity being the other. Using complete Galerkin projection 
with Pi = (2,2,2),i = 1,2, gives rise to a 60 state ODE model. To apply 
WR/AWR to this system, we first identify the weakly interacting states. 
By construction the two zones weakly affect each other, which is identified 
by the spectral clustering |13j (or wave equation based clustering [1^) ap- 



plied to the system (94). This decomposition is imposed on the complete 
Galerkin system, as explained in section 4.3.1 As expected, we found that 



if ones applies spectral clustering to the complete Galerkin system instead, 
one recovers same decomposition. 

Treating 1000 Monte Carlo samples as the truth, we compare the results 
of simulated full Galerkin projected system using both standard waveform 
relaxation [TH] as well as adaptive waveform relaxation |17] in Fig. [5^). AWR 
provides a speed-up by a factor of ^ 12. In Fig. [5^), one can visually see that 
the complete Galerkin Projection predicts the same temperature variation 
over 8 hours as Monte Carlo based methods. 

As explained before, one can further exploit the weak interaction between 
the two zones to reduce the overall number of equations in Galerkin projec- 
tion. To construct the AGP system, we reduce the order of expansion for 
the random parameters indirectly affecting each zone so that Pi = (2,2, 1) 
and P2 = (1, 2, 2). With this the number of equations in Galerkin projection 
reduces from 60 to 50. The resulting solution is shown in Fig. ^). We see 
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that the error starts to grow as time increases. However, over 8 hours the 
max error in the room temperatures is 5 x lO^^K. Thus, despite reducing 
the computational effort, one can still get a fairly accurate answer. 



Figure [Sp) shows the effect of coupling (which is the reciprocal of the 

coefficient of thermal conductivity of the internal wall) on errors introduced 
in complete and approximate Galerkin projections. As expected, the ap- 
proximate Galerkin projection has higher error (given by i?T(i)) than com- 
plete Galerkin projection (given by Ec{t)). Moreover this error is more 
pronounced at low iteration numbers. From the figure, it also clear that 
as the coupling increases, the number of iterations required for obtaining 
same solution accuracy increases. For further discussion on the relationship 
between the coupling and number of iterations, see |17j . 

5.2.2 Multi Zone Example 

In this section, we consider a larger 6 zone building thermal network model 
with 68 states. This model admits a decomposition into 23 subsystems, as 
revealed by the spectral graph approach (see figure [6}o) . This decomposition 
is consistent with three different time scales (associated with external and 
internal wall temperature, and internal zone temperatures) present in the 
system, as shown by the three bands in figure^). 




Next we demonstrate non-intrusive PWR approach to compute uncer- 
tainty in energy consumption due to both parametric uncertainty and time 
varying uncertain loads. As described earlier, we use KL expansion to trans- 
form time varying uncertainty into parametric form. 

KL Expansion [25j: Let {Xt = X{S^,t),t E [a, 6]} be a quadratic mean 
square second-order stochastic process with covariance functions R{t,s). If 
{(pn{t)} are the eigenfunctions of the integral operator with kernel R{-, •) 
and {A„} the corresponding eigenvalues, i.e. 






t e [a, b] 



(95) 



then. 



N 



X{t,0) = X{t) + lim \/^an{C)4>n{t), uniformaly for t G [a, 6] 



n=l 



(96) 
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where, X{t) is the mean of the process and the hmit is taken in the quadratic 
mean sense. The random coefficients {a„} satisfy 



1 



'Xr, 



b 

iXiC,t)-X{t))Mt)dt (97) 



and are uncorrelated E[aman] = ^mn- The basis functions also satisfy or- 
thogonahty property 

fb 

(t>m{t)Mt)dt = 5mn. (98) 



and the kernel admits an expansion of the form 

N 

Af-s>oo 



R{s,t)= lim VA„(/>„(t)0n(s). (99) 

N — ^ — * 



n=l 



Generally, analytical solution to the eigenvalue problem (95), also known 
as Fredholm equation of second kind is not available. Several numerical 
techniques have been proposed, we used the expansion method described in 

m- 

For applying KL expansion to the building problem, we assume that 
the stochastic disturbances {Tamb{i),Qs{t)-,Qint{t)) are Gaussian processes. 
This guarantees that the random variables a„ in the KL expansion are inde- 
pendent Gaussian random variables with a zero mean ([26]). Let the joint 
distribution of a nonstationary Gaussian process be, 

fiXi X ) = - e 2{l-p2(t,s)) \<T-^(ty a^is) <T{t)<T{B) 

' 27ra{s)a{t)^l-p{t,s) 

(100) 

where, p{t, s) is the correlation coefficient and is related to covariance kernel 
as 

R{t, s) = p{t, s)a{t)a{s). (101) 

We assumed the processes Tamb{t),Qs{t) to have a stationary exponential 
correlation function 

R{t,s) = a'^e"^ , (102) 

with a constant variance cr^ and a constant correlation time scale Tg. For 
the internal occupancy load Qint{t) we constructed i?(t, s) as follows. For a 
typical office building, we know that the occupancy load is negligible with 
low variance during early and later parts of the day. On the other hand 
during peak hours in the middle of the day the occupant load can show 
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significantly high variabihty. To capture this effect we divided the normal- 
ized time domain [0, 1] = [0,ti] U (ti,t2) U (t2, 1) and obtained the desired 



variation by choosing (in expression 101 ) 



\t-s\ 

fT(t) = cr(tanh(a(t - ti)) - tanh(a(x - t2)))/2, s) = e Tc(s,t) ^ (103) 

with the correlation time scale 

Tc{s,t) = (l-tanh(a(t-ti)))(l-tanh(a(s-ti)))/4 

+ (1 + tanh(a(t - t2)))(l + tanh(a(s - t2)))/4, (104) 

and parameter a controls the slope of tanh function. Figure [7] shows the 
covariance kernel for external (Tamb{'t),Qsit)) and internal Qint{t) loads. 
For the choice of parameters indicated in the figure [7j we found using the 
expansion method p6] with Legendre polynomials as the basis functions, 
that KL expansion upto order 3 and upto order 6 can capture more that 
90% of total variance, for internal and external loads, respectively. 

In UQ computation, we considered the effect of 14 random variables com- 
prising of external wall thermal resistance in the 6 zones, and first dominant 
random variable obtained in the KL representation of internal load (for each 
zone) and first two dominant random variable obtained in KL expansion for 
solar load. Figure [8] show the non-intrusive PWR results on the decomposed 
network model. As is evident, the iterations converge rapidly in two steps 
with a distribution close to that obtained from QMC (using a 25000-sample 



Sobol sequence) applied to the 68 thermal network model (94) as a whole. 



5.3 Coupled Oscillators 

Finally, we consider a coupled phase only oscillator system which is governed 
by nonlinear equations 

N 

Xi = u}i + '^Kijsm{xj - Xi), z = l,---,n, (105) 

i=i 

where, n = 80 is the number of oscillators, cjj, i = 1, • • • , n is the angular 
frequency of oscillators and K = [Kij] is the coupling matrix. The frequen- 
cies uji of every alternative oscillator i.e. i = 1, 3, • • • , 79 is assumed to be 
uncertain with a Gaussian distribution with 20% tolerance (i.e. with a total 
p = 40 uncertain parameters); all the other parameters are assumed to take 
a fixed mean value. We are interested in the distribution of the synchroniza- 
tion parameters, R{t) and phase (/)(t), defined by i?(t)e'^'^*) = Ylf=i e^^^^^\ 
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Figure [o] shows the topology of the network of oscihators (left panel) , along 
with the eigenvalue spectrum of the graph Laplacian (right panel). The 
spectral gap at 40, implies 40 weakly interacting subsystems in the network. 

Figure 10 shows UQ results obtained by application of non-intrusive 
PWR to the decomposed system with Ig = 5, Ic = 2. We make comparison 
with QMC, in which the complete system (105) is solved at 25,000 Sobol 
points [27J. Remarkably the PWR converges in 4 — 5 iterations giving very 
similar results to that of QMC. It would be infeasible to use full grid colloca- 
tion for the networks as a whole, since even with lowest level of collocation 
grid, i.e. / = 2 for each parameter, the number of samples required become 

= 2^0 = 1.0995e + 012!. 



6 Conclusion and Future Work 

In this paper we have proposed an uncertainty quantification approach which 
exploits the underlying dynamics and structure of the system. In specific we 
considered a class of networked system, whose subsystems are dynamically 
weakly coupled to each other. We showed how these weak interactions can be 
exploited to overcome the dimensionality curse associated with traditional 
UQ methods. By integrating graph decomposition and waveform relaxation 
with generalized polynomial chaos and probabilistic collocation framework, 
we proposed an iterative UQ approach which we called probabilistic wave- 
form relaxation. We developed both intrusive and non-intrusive forms of 
PWR. We proved that this iterative scheme converges and illustrated it on 
several examples with promising results. Several questions need to be fur- 
ther investigated, these include: how the choice parameters associated with 
PWR algorithm affects its rate of convergence and the approximation error. 
In order to exploit multiple time scales that may be present in a system, 
multigrid extension [28j of PWR will be desirable. 
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Figure 5: (a) Comparison of Monte Carlo, complete Galerkin projection and 
approximate Galerkin projection. (b)Normalized error in waveform relax- 
ation as a function of iteration count with increasing coupling. Complete 
Galerkin and approximate Galerkin are shown. Approximate Galerkin sys- 
tem is found to have greater error as a function of iteration number. 
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Figure 6: (a) Shows the there bands of eigenvalues of the time averaged 
A(t; ^) for nominal parameter values, (b) First spectral gap in graph Lapla- 
cian revealing 23 subsystems in the network model. 
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Iteralion 1 , PWR Iteration 2, PWR QMC 




kW-hr kW-hr kW-hr 

Figure 8: Histogram of building energy computation for two iterations in 
PWR. Also shown is the corresponding histogram obtained by QMC for 
comparison. 




Figure 9: Left panel shows a network of iV = 80 phase only oscillators. Right 
panel shows spectral gap in eigenvalues of normalized graph Laplacian, that 
reveals that there are 40 weakly interacting subsystems. 
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Figure 10: Convergence of mean of the magnitude R{t) and phase (j){t), and 
the respective histograms at t = 0.5. 
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